
include("../src/qmcmary.jl")
using ..qmcmary

using Test
using ReverseDiff: jacobian, jacobian!
using LinearAlgebra

BLAS.set_num_threads(1)


"""
不使用稳定
"""
function rawgbar(allbmats, idx)
    Nx = size(allbmats[1])[1]
    Rs = I(Nx)
    for i in Base.OneTo(idx)
        Rs = allbmats[i]*Rs
    end
    Ls = Diagonal(ones(Nx))
    for i in length(allbmats):-1:(idx+1)
        Ls = Ls*allbmats[i]
    end
    return inv(Diagonal(ones(Nx)) + inv(Ls)*inv(Rs))
end


function thetabar_example()
    #通过一个[t1, ..., tn]在确定的hscfg下计算出sign和 d ln s / dti
    #微调ti计算sign，验证数值正确性
    Nx = 2
    Nt = 20
    Ng = 5
    sslen = Int(Nt//Ng)
    bmats1 = Vector{Matrix{ComplexF64}}(undef, Ng)
    bmats2 = Vector{Matrix{ComplexF64}}(undef, Ng)
    allbmats1 = Vector{Matrix{ComplexF64}}(undef, Nt)
    allbmats2 = Vector{Matrix{ComplexF64}}(undef, Nt)
    bidxs = Vector{Int}(undef, Ng)
    hscfg = Matrix{Int}(undef, Nt, Nx)
    psi1 = rand(Nx)
    th1 = rand(Nx)
    # ∂ B_tau_xy  / ∂ nz，实部和虚部
    pbrpnx = Array{Float64, 4}(undef, Nt, 2*Nx, 2*Nx, Nx)
    pbipnx = Array{Float64, 4}(undef, Nt, 2*Nx, 2*Nx, Nx)
    pbrpny = Array{Float64, 4}(undef, Nt, 2*Nx, 2*Nx, Nx)
    pbipny = Array{Float64, 4}(undef, Nt, 2*Nx, 2*Nx, Nx)
    pbrpnz = Array{Float64, 4}(undef, Nt, 2*Nx, 2*Nx, Nx)
    pbipnz = Array{Float64, 4}(undef, Nt, 2*Nx, 2*Nx, Nx)
    #
    nx1 = sin.(psi1).*sin.(th1)
    ny1 = cos.(psi1).*sin.(th1)
    nz1 = cos.(th1)
    nx3 = copy(nx1)
    ny3 = copy(ny1)
    nz3 = copy(nz1)
    ny3[2] += 0.01 
    #
    hk = lattice_chain(Nx, -1.0)#[0.0 1.0; 1.0 0.0]
    pet = rand(ComplexF64, Nx, Nx)
    pet = pet + adjoint(pet)
    hk += pet
    Ui = 10*rand(Nx)
    dt = 0.1
    ehk = exp(-dt*hk)
    ehk = kron([1 0; 0 1], ehk)
    for bi in Base.OneTo(Ng)
        cfg = 2(round.(rand(Nx)) .- 0.5)
        hscfg[bi, :] .= cfg
        #println(bi)
        bidxs[bi] = bi
        bmr, bmi = bmat_IsingAD(ehk, cfg, Ui, nx1, ny1, nz1, dt)
        #println(bmr)
        bmats1[bi] = complex.(bmr, bmi)
        allbmats1[bi] = bmats1[bi]
        bmr, bmi = bmat_IsingAD(ehk, cfg, Ui, nx3, ny3, nz3, dt)
        #println(bmr)
        bmats2[bi] = complex.(bmr, bmi)
        allbmats2[bi] = bmats2[bi]
        #求 ∂ B_tau_xy  / ∂ nx1, ∂ B_tau_xy  / ∂ ny1
        _fbr(nx, ny, nz) = bmat_IsingAD(ehk, cfg, Ui, nx, ny, nz, dt)[1]
        _fbi(nx, ny, nz) = bmat_IsingAD(ehk, cfg, Ui, nx, ny, nz, dt)[2]
        pbpnxr, pbpnyr, pbpnzr = jacobian(_fbr, (nx1, ny1, nz1))
        pbrpnx[bi, :, :, :] = reshape(pbpnxr, 2*Nx, 2*Nx, Nx)
        pbrpny[bi, :, :, :] = reshape(pbpnyr, 2*Nx, 2*Nx, Nx)
        pbrpnz[bi, :, :, :] = reshape(pbpnzr, 2*Nx, 2*Nx, Nx)
        pbpnxi, pbpnyi, pbpnzi = jacobian(_fbi, (nx1, ny1, nz1))
        pbipnx[bi, :, :, :] = reshape(pbpnxi, 2*Nx, 2*Nx, Nx)
        pbipny[bi, :, :, :] = reshape(pbpnyi, 2*Nx, 2*Nx, Nx)
        pbipnz[bi, :, :, :] = reshape(pbpnzi, 2*Nx, 2*Nx, Nx)
        #println(pbrpnz[bi, :, :, 1])
        #println(pbipnz[bi, :, :, 1])
        #exit()
    end
    ss1 = ScrollSVD(bmats1)
    ss2 = ScrollSVD(bmats2)
    for si in Base.OneTo(sslen-1)
        for bi in Base.OneTo(Ng)
            cfg = 2(round.(rand(Nx)) .- 0.5)
            hscfg[bi+si*Ng, :] .= cfg
            #println(bi+si*Ng)
            bidxs[bi] = bi+si*Ng
            bmr, bmi = bmat_IsingAD(ehk, cfg, Ui, nx1, ny1, nz1, dt)
            #println(bmi)
            bmats1[bi] = complex.(bmr, bmi)
            allbmats1[bi+si*Ng] = bmats1[bi]
            bmr, bmi = bmat_IsingAD(ehk, cfg, Ui, nx3, ny3, nz3, dt)
            #println(bmi)
            bmats2[bi] = complex.(bmr, bmi)
            allbmats2[bi+si*Ng] = bmats2[bi]
            #
            _fbr(nx, ny, nz) = bmat_IsingAD(ehk, cfg, Ui, nx, ny, nz, dt)[1]
            _fbi(nx, ny, nz) = bmat_IsingAD(ehk, cfg, Ui, nx, ny, nz, dt)[2]
            #
            pbpnxr, pbpnyr, pbpnzr = jacobian(_fbr, (nx1, ny1, nz1))
            pbrpnx[bi+si*Ng, :, :, :] = reshape(pbpnxr, 2*Nx, 2*Nx, Nx)
            pbrpny[bi+si*Ng, :, :, :] = reshape(pbpnyr, 2*Nx, 2*Nx, Nx)
            pbrpnz[bi+si*Ng, :, :, :] = reshape(pbpnzr, 2*Nx, 2*Nx, Nx)
            #println(pbrpnz[bi+si*Ng, :, :, :])
            pbpnxi, pbpnyi, pbpnzi = jacobian(_fbi, (nx1, ny1, nz1))
            pbipnx[bi+si*Ng, :, :, :] = reshape(pbpnxi, 2*Nx, 2*Nx, Nx)
            pbipny[bi+si*Ng, :, :, :] = reshape(pbpnyi, 2*Nx, 2*Nx, Nx)
            pbipnz[bi+si*Ng, :, :, :] = reshape(pbpnzi, 2*Nx, 2*Nx, Nx)
            #println(pbipny[bi+si*Ng, :, :, :])
        end
        push!(ss1, bmats1)
        push!(ss2, bmats2)
    end
    #
    #for bi in Base.OneTo(Nt)
    #    if hscfg[bi, 1] == 1 && hscfg[bi, 2] == 1
    #        println(pbipny[bi, 1, :, 1], pbipny[bi, 2, :, 2])
    #        println(pbipny[bi, 3, :, 1], pbipny[bi, 4, :, 2])
    #        println(pbrpnz[bi, 1, :, 1], pbrpnz[bi, 2, :, 2])
    #        println(pbrpnz[bi, 3, :, 1], pbrpnz[bi, 4, :, 2])
    #        @assert false
    #    end
    #end
    #
    sp = default_splitting(Nt, hk, Ui)
    tapemap = pbpn_calc_tapemap(sp, nx1, ny1, nz1)
    pbrpnx2, pbipnx2, pbrpny2, pbipny2, pbrpnz2, pbipnz2 = pbpn_calc(sp, hscfg, nx1, ny1, nz1; tapemap=tapemap)
    println(size(pbrpny), size(pbrpny2))
    println(size(isapprox.(pbrpny, pbrpny2, atol=1e-13)))
    @testset "∂B/∂n" begin
        @test all(isapprox.(pbrpnx, pbrpnx2, atol=1e-14))
        @test all(isapprox.(pbrpny, pbrpny2, atol=1e-14))
        @test all(isapprox.(pbrpnz, pbrpnz2, atol=1e-14))
        @test all(isapprox.(pbipnx, pbipnx2, atol=1e-14))
        @test all(isapprox.(pbipny, pbipny2, atol=1e-14))
        @test all(isapprox.(pbipnz, pbipnz2, atol=1e-14))
    end
    # 求 \bar  B_tau_xy
    bbar2 = Array{ComplexF64, 3}(undef, Nt, 2*Nx, 2*Nx)
    for bi in Base.OneTo(Nt)
        gs1 = rawgbar(allbmats1, bi-1)
        bba = bbar_from_gbar(gs1, allbmats1[bi])
        bbar2[bi, :, :] = bba
    end
    bbar = bbar_calc(ss1, allbmats1)
    @test all(isapprox.(bbar2, bbar, atol=1e-8))
    # \bar nx1 = ∑_{tau,xy} \bar  B_tau_xy ∂ B_tau_xy  / ∂ nx1
    nxbar = zeros(Nx)
    nybar = zeros(Nx)#Vector{Float64}(undef, Nx)
    nzbar = zeros(Nx)#Vector{Float64}(undef, Nx)
    #for si in Base.OneTo(Nx)
    #注意虚部附加的负号
    for bi in Base.OneTo(Nt); for xi in Base.OneTo(2*Nx); for yi in Base.OneTo(2*Nx)
        nzbar += pbrpnz[bi, xi, yi, :] * real(bbar[bi, xi, yi])
        nzbar -= pbipnz[bi, xi, yi, :] * imag(bbar[bi, xi, yi])
        nybar += pbrpny[bi, xi, yi, :] * real(bbar[bi, xi, yi])
        nybar -= pbipny[bi, xi, yi, :] * imag(bbar[bi, xi, yi])
        nxbar += pbrpnx[bi, xi, yi, :] * real(bbar[bi, xi, yi])
        nxbar -= pbipnx[bi, xi, yi, :] * imag(bbar[bi, xi, yi])
    end; end; end
    #end
    #println(nxbar, nybar, nzbar)
    tapemap = pbpn_calc_tapemap(sp, nx1, ny1, nz1)
    nxbar2, nybar2, nzbar2 = meas_grad(ss1, allbmats1, sp, hscfg, nx1, ny1, nz1; tapemap=tapemap)
    @testset "n̄" begin
        @test all(isapprox.(nxbar, nxbar2, atol=1e-13))
        @test all(isapprox.(nybar, nybar2, atol=1e-13))
        @test all(isapprox.(nzbar, nzbar2, atol=1e-13))
    end
    #求sign
    sig1 = real(det(I(2*Nx) + prod(reverse(allbmats1))))
    println(sig1)
    println(-log(sig1))
    sig2 = real(det(I(2*Nx) + prod(reverse(allbmats2))))
    println(sig2)
    println(-log(sig2))
    println(-log(sig2) + log(sig1))
    println(nybar)
    #求∂ny/∂θ和∂nz/∂θ
    pnxpt = sin.(psi1) .* cos.(th1)
    pnypt = cos.(psi1) .* cos.(th1)
    pnzpt = -sin.(th1)
    pnxpp = cos.(psi1) .* sin.(th1)
    pnypp = -sin.(psi1) .* sin.(th1)
    pnzpp = zeros(Nx)
    #求∂L/∂θ
    thetabar = @. nxbar * pnxpt + nybar * pnypt +  nzbar * pnzpt
    println(thetabar)
    psibar = @. nxbar * pnxpp + nybar * pnypp +  nzbar * pnzpp
    println(psibar)
end



function sgn_scratch(ss::ScrollSVD{T}) where T
    siz = size(ss.B[end])
    ptr = findfirst(ss.L)
    if isnothing(ptr)
        VL = Diagonal(ones(siz[1]))
        DL = VL
        UL = VL
        UR, DR, VR = ss.F[end].U, Diagonal(ss.F[end].S), ss.F[end].Vt
    elseif ptr == 1
        VL, DL, UL = ss.F[1].U, Diagonal(ss.F[1].S), ss.F[1].Vt
        UR = Diagonal(ones(siz[1]))
        DR = UR
        VR = UR
        ##上面这个数值稳定性会突然出问题
        ##用下面这个
        ##VL DL UL UR=I DR=I VR=I -> I I I UR=VL DR=DL VR=UL
        #VL = Diagonal(ones(siz[1]))
        #DL = VL
        #UL = VL
        #UR, DR, VR = ss.F[1].U, Diagonal(ss.F[1].S), ss.F[1].Vt
    else
        VL, DL, UL = ss.F[ptr].U, Diagonal(ss.F[ptr].S), ss.F[ptr].Vt
        UR, DR, VR = ss.F[ptr-1].U, Diagonal(ss.F[ptr-1].S), ss.F[ptr-1].Vt
    end
    #gtt = inv(Diagonal(ones(siz[1]))+UR*DR*VR*VL*DL*UL)
    #M = inv(UL*UR) + DR*(VR*VL)*DL
    #Fm = svd(M)
    #gtt = inv(Fm.Vt*UL)*inv(Diagonal(Fm.S))*inv(UR*Fm.U)
    #
    DLS = Diagonal(ones(Float64, siz[1]))
    DLB = Diagonal(ones(Float64, siz[1]))
    DRS = Diagonal(ones(Float64, siz[1]))
    DRB = Diagonal(ones(Float64, siz[1]))
    for i in Base.OneTo(siz[1])
        if DL[i, i] > 1.0
            DLB[i, i] = DL[i, i]
            DLS[i, i] = 1.0
        else
            DLS[i, i] = DL[i, i]
            DLB[i, i] = 1.0
        end
        if DR[i, i] > 1.0
            DRB[i, i] = DR[i, i]
            DRS[i, i] = 1.0
        else
            DRS[i, i] = DR[i, i]
            DRB[i, i] = 1.0
        end
    end
    #
    M = inv(DRB)*adjoint(UL*UR)*inv(DLB) + DRS*(VR*VL)*DLS
    Fm = svd(M, alg=LinearAlgebra.QRIteration())
    #ML = adjoint(UL)*inv(DLB)*adjoint(Fm.Vt)
    #MR = adjoint(Fm.U)*inv(DRB)*adjoint(UR)
    #gtt = ML*inv(Diagonal(Fm.S))*MR
    #增加计算phase
    sgn = det(Fm.Vt*UL)*det(UR*Fm.U)
    sgn = log(abs(sgn))
    for sval in Fm.S
        sgn += log(sval)
    end
    for sval in diag(DLB)
        sgn += log(sval)
    end
    for sval in diag(DRB)
        sgn += log(sval)
    end
    return sgn
end

function thetabar_example2()
    #通过一个[t1, ..., tn]在确定的hscfg下计算出sign和 d ln s / dti
    #微调ti计算sign，验证数值正确性
    L = 3
    Nx = 3*L^2
    psi = 1.7652
    the = -1.0
    nx = sin(psi)*sin(the)*ones(Nx)
    ny = cos(psi)*sin(the)*ones(Nx)
    nz = cos(the)*ones(Nx)
    #
    hk = lattice_kagome(ComplexF64, L, -1.0+0.0im)
    hk = kron([1 0; 0 1], hk)
    #lattice_hexagonal(ComplexF64, L, -1.0+0.0im; λ=1.0+0.0im)
    pet = rand(ComplexF64, 2*Nx, 2*Nx)
    pet = pet + adjoint(pet)
    hk += pet
    Ui = 6*ones(Nx)
    Nt = 50
    Ng = 5
    hscfg = Matrix{Int}(undef, Nt, Nx)
    for bi in Base.OneTo(Nt)
        cfg = 2(round.(rand(Nx)) .- 0.5)
        hscfg[bi, :] .= cfg
    end
    sp = default_splitting(Nt, hk, Ui; Z2=false)
    sslen, bmats, allbmats, ss = initialize_SS(Nt, Ng, sp, hscfg, nx, ny, nz)
    logS1 = sgn_scratch(ss)
    tapemap = pbpn_calc_tapemap(sp, nx, ny, nz)
    nxbar, nybar, nzbar = meas_grad(ss, allbmats, sp, hscfg, nx, ny, nz; tapemap=tapemap)
    #println(nybar)
    #println(nzbar)
    #
    admat = zeros(Float64, Nx, 5)
    ndmat = zeros(Float64, Nx, 5)
    for xi in Base.OneTo(Nx)
        nx1 = copy(nx)
        ny1 = copy(ny)
        nz1 = copy(nz)
        nx1[xi] += 0.001
        sslen, bmats, allbmats, ss = initialize_SS(Nt, Ng, sp, hscfg, nx1, ny1, nz1)
        logS2 = sgn_scratch(ss)
        #d log S = - d logabs(w)
        #所以这里需要加一个负号
        #println(-logS2 + logS1)
        nd = -logS2 + logS1
        ad = 0.001*nxbar[xi]
        #println(nd, " ", ad)
        admat[xi, 1] = ad
        ndmat[xi, 1] = nd
        #@test isapprox(nd, ad, rtol=1e-2, atol=2e-4)
        #println(0.001*nxbar[10])
        nx1[xi] -= 0.001
        ny1[xi] += 0.001
        sslen, bmats, allbmats, ss = initialize_SS(Nt, Ng, sp, hscfg, nx1, ny1, nz1)
        logS2 = sgn_scratch(ss)
        #d log S = - d logabs(w)
        #所以这里需要加一个负号
        #println(-logS2 + logS1)
        nd = -logS2 + logS1
        ad = 0.001*nybar[xi]
        #println(nd, " ", ad)
        admat[xi, 2] = ad
        ndmat[xi, 2] = nd
        #@test isapprox(nd, ad, rtol=1e-2, atol=2e-4)
        #
        ny1[xi] -= 0.001
        nz1[xi] += 0.001
        sslen, bmats, allbmats, ss = initialize_SS(Nt, Ng, sp, hscfg, nx1, ny1, nz1)
        logS2 = sgn_scratch(ss)
        #d log S = - d logabs(w)
        #所以这里需要加一个负号
        #println(-logS2 + logS1)
        nd = -logS2 + logS1
        ad = 0.001*nzbar[xi]
        #println(nd, " ", ad)
        admat[xi, 3] = ad
        ndmat[xi, 3] = nd
        #@test isapprox(nd, ad, rtol=1e-2, atol=2e-4)
        #
        pnxpt = sin.(psi) .* cos.(the)
        pnypt = cos.(psi) .* cos.(the)
        pnzpt = -sin.(the)
        pnxpp = cos.(psi) .* sin.(the)
        pnypp = -sin.(psi) .* sin.(the)
        pnzpp = zeros(Nx)
        #
        psi1 = psi*ones(Nx)
        the1 = the*ones(Nx)
        psi1[xi] += 0.001
        nx1 = @. sin(psi1)*sin(the1)
        ny1 = @. cos(psi1)*sin(the1)
        nz1 = @. cos(the1)
        sslen, bmats, allbmats, ss = initialize_SS(Nt, Ng, sp, hscfg, nx1, ny1, nz1)
        logS2 = sgn_scratch(ss)
        psibar = @. nxbar * pnxpp + nybar * pnypp +  nzbar * pnzpp
        nd = -logS2 + logS1
        ad = 0.001*psibar[xi]
        #println(nd, " ", ad)
        admat[xi, 4] = ad
        ndmat[xi, 4] = nd
        #@test isapprox.(nd, ad, rtol=1e-2, atol=2e-4)
        #
        psi1 = psi*ones(Nx)
        the1 = the*ones(Nx)
        the1[xi] += 0.001
        nx1 = @. sin(psi1)*sin(the1)
        ny1 = @. cos(psi1)*sin(the1)
        nz1 = @. cos(the1)
        sslen, bmats, allbmats, ss = initialize_SS(Nt, Ng, sp, hscfg, nx1, ny1, nz1)
        logS2 = sgn_scratch(ss)
        thetabar = @. nxbar * pnxpt + nybar * pnypt +  nzbar * pnzpt
        nd = -logS2 + logS1
        ad = 0.001*thetabar[xi]
        #println(nd, " ", ad)
        admat[xi, 5] = ad
        ndmat[xi, 5] = nd
        #@test isapprox.(nd, ad, rtol=1e-2, atol=2e-4)
    end
    @testset "-∂log(s)/∂x" begin
        @test all(isapprox.(ndmat, admat, rtol=1e-2, atol=2e-4))
    end
end


function thetabar_example3()
    #通过一个[t1, ..., tn]在确定的hscfg下计算出sign和 d ln s / dti
    #微调ti计算sign，验证数值正确性
    Nx = 6
    Nt = 40
    Ng = 5
    sslen = Int(Nt//Ng)
    bmats1 = Vector{Matrix{ComplexF64}}(undef, Ng)
    bmats2 = Vector{Matrix{ComplexF64}}(undef, Ng)
    allbmats1 = Vector{Matrix{ComplexF64}}(undef, Nt)
    allbmats2 = Vector{Matrix{ComplexF64}}(undef, Nt)
    bidxs = Vector{Int}(undef, Ng)
    hscfg = Matrix{Int}(undef, Nt, Nx)
    psi1 = rand(Nx)
    th1 = rand(Nx)
    # ∂ B_tau_xy  / ∂ nz，实部和虚部
    pbrpnx = Array{Float64, 4}(undef, Nt, 2*Nx, 2*Nx, Nx)
    pbipnx = Array{Float64, 4}(undef, Nt, 2*Nx, 2*Nx, Nx)
    pbrpny = Array{Float64, 4}(undef, Nt, 2*Nx, 2*Nx, Nx)
    pbipny = Array{Float64, 4}(undef, Nt, 2*Nx, 2*Nx, Nx)
    pbrpnz = Array{Float64, 4}(undef, Nt, 2*Nx, 2*Nx, Nx)
    pbipnz = Array{Float64, 4}(undef, Nt, 2*Nx, 2*Nx, Nx)
    #
    nx1 = sin.(psi1).*sin.(th1)
    ny1 = cos.(psi1).*sin.(th1)
    nz1 = cos.(th1)
    nx3 = copy(nx1)
    ny3 = copy(ny1)
    nz3 = copy(nz1)
    ny3[2] += 0.01 
    #
    hk = lattice_chain(Nx, -1.0)#[0.0 1.0; 1.0 0.0]
    pet = rand(ComplexF64, Nx, Nx)
    pet = pet + adjoint(pet)
    hk += pet
    Ui_ = 3*rand(Nx)
    sp = default_splitting3(Nt, hk, Ui_)
    dt = 0.1
    ehk1 = exp(-(dt+dt*im)*hk)
    ehk1 = kron([1 0; 0 1], ehk1)
    ehk2 = exp(-(dt-dt*im)*hk)
    ehk2 = kron([1 0; 0 1], ehk2)
    for bi in Base.OneTo(Ng)
        cfg = 2(round.(rand(Nx)) .- 0.5)
        hscfg[bi, :] .= cfg
        #println(bi)
        bidxs[bi] = bi
        ehk, Ui = unpack_splitting(sp, bi)
        if bi % 2 == 1
            @assert all(isapprox.(ehk, ehk1, atol=1e-14))
            @assert all(isapprox.(Ui, 2*Ui_, atol=1e-14))
        else
            @assert all(isapprox.(ehk, ehk2, atol=1e-14))
            @assert all(isapprox.(Ui, zeros(Nx), atol=1e-14))
        end
        bmr, bmi = bmat_IsingAD(ehk, cfg, Ui, nx1, ny1, nz1, dt)
        #println(bmr)
        bmats1[bi] = complex.(bmr, bmi)
        allbmats1[bi] = bmats1[bi]
        bmr, bmi = bmat_IsingAD(ehk, cfg, Ui, nx3, ny3, nz3, dt)
        #println(bmr)
        bmats2[bi] = complex.(bmr, bmi)
        allbmats2[bi] = bmats2[bi]
        #求 ∂ B_tau_xy  / ∂ nx1, ∂ B_tau_xy  / ∂ ny1
        _fbr(nx, ny, nz) = bmat_IsingAD(ehk, cfg, Ui, nx, ny, nz, dt)[1]
        _fbi(nx, ny, nz) = bmat_IsingAD(ehk, cfg, Ui, nx, ny, nz, dt)[2]
        pbpnxr, pbpnyr, pbpnzr = jacobian(_fbr, (nx1, ny1, nz1))
        pbrpnx[bi, :, :, :] = reshape(pbpnxr, 2*Nx, 2*Nx, Nx)
        pbrpny[bi, :, :, :] = reshape(pbpnyr, 2*Nx, 2*Nx, Nx)
        pbrpnz[bi, :, :, :] = reshape(pbpnzr, 2*Nx, 2*Nx, Nx)
        pbpnxi, pbpnyi, pbpnzi = jacobian(_fbi, (nx1, ny1, nz1))
        pbipnx[bi, :, :, :] = reshape(pbpnxi, 2*Nx, 2*Nx, Nx)
        pbipny[bi, :, :, :] = reshape(pbpnyi, 2*Nx, 2*Nx, Nx)
        pbipnz[bi, :, :, :] = reshape(pbpnzi, 2*Nx, 2*Nx, Nx)
        #println(pbrpnz[bi, :, :, 1])
        #println(pbipnz[bi, :, :, 1])
        #exit()
    end
    ss1 = ScrollSVD(bmats1)
    ss2 = ScrollSVD(bmats2)
    for si in Base.OneTo(sslen-1)
        for bi in Base.OneTo(Ng)
            cfg = 2(round.(rand(Nx)) .- 0.5)
            hscfg[bi+si*Ng, :] .= cfg
            #println(bi+si*Ng)
            bidxs[bi] = bi+si*Ng
            ehk, Ui = unpack_splitting(sp, bi+si*Ng)
            if (bi+si*Ng) % 2 == 1
                @assert all(isapprox.(ehk, ehk1, atol=1e-14))
            else
                @assert all(isapprox.(ehk, ehk2, atol=1e-14))
            end
            bmr, bmi = bmat_IsingAD(ehk, cfg, Ui, nx1, ny1, nz1, dt)
            #println(bmi)
            bmats1[bi] = complex.(bmr, bmi)
            allbmats1[bi+si*Ng] = bmats1[bi]
            bmr, bmi = bmat_IsingAD(ehk, cfg, Ui, nx3, ny3, nz3, dt)
            #println(bmi)
            bmats2[bi] = complex.(bmr, bmi)
            allbmats2[bi+si*Ng] = bmats2[bi]
            #
            _fbr(nx, ny, nz) = bmat_IsingAD(ehk, cfg, Ui, nx, ny, nz, dt)[1]
            _fbi(nx, ny, nz) = bmat_IsingAD(ehk, cfg, Ui, nx, ny, nz, dt)[2]
            #
            pbpnxr, pbpnyr, pbpnzr = jacobian(_fbr, (nx1, ny1, nz1))
            pbrpnx[bi+si*Ng, :, :, :] = reshape(pbpnxr, 2*Nx, 2*Nx, Nx)
            pbrpny[bi+si*Ng, :, :, :] = reshape(pbpnyr, 2*Nx, 2*Nx, Nx)
            pbrpnz[bi+si*Ng, :, :, :] = reshape(pbpnzr, 2*Nx, 2*Nx, Nx)
            #println(pbrpnz[bi+si*Ng, :, :, :])
            pbpnxi, pbpnyi, pbpnzi = jacobian(_fbi, (nx1, ny1, nz1))
            pbipnx[bi+si*Ng, :, :, :] = reshape(pbpnxi, 2*Nx, 2*Nx, Nx)
            pbipny[bi+si*Ng, :, :, :] = reshape(pbpnyi, 2*Nx, 2*Nx, Nx)
            pbipnz[bi+si*Ng, :, :, :] = reshape(pbpnzi, 2*Nx, 2*Nx, Nx)
            #println(pbipny[bi+si*Ng, :, :, :])
        end
        push!(ss1, bmats1)
        push!(ss2, bmats2)
    end
    #
    #for bi in Base.OneTo(Nt)
    #    if hscfg[bi, 1] == 1 && hscfg[bi, 2] == 1
    #        println(pbipny[bi, 1, :, 1], pbipny[bi, 2, :, 2])
    #        println(pbipny[bi, 3, :, 1], pbipny[bi, 4, :, 2])
    #        println(pbrpnz[bi, 1, :, 1], pbrpnz[bi, 2, :, 2])
    #        println(pbrpnz[bi, 3, :, 1], pbrpnz[bi, 4, :, 2])
    #        @assert false
    #    end
    #end
    #
    tapemap = pbpn_calc_tapemap(sp, nx1, ny1, nz1)
    pbrpnx2, pbipnx2, pbrpny2, pbipny2, pbrpnz2, pbipnz2 = pbpn_calc(sp, hscfg, nx1, ny1, nz1; tapemap=tapemap)
    println(size(pbrpny), size(pbrpny2))
    println(size(isapprox.(pbrpny, pbrpny2, atol=1e-13)))
    @testset "∂B/∂n" begin
        @test all(isapprox.(pbrpnx, pbrpnx2, atol=1e-14))
        @test all(isapprox.(pbrpny, pbrpny2, atol=1e-14))
        @test all(isapprox.(pbrpnz, pbrpnz2, atol=1e-14))
        @test all(isapprox.(pbipnx, pbipnx2, atol=1e-14))
        @test all(isapprox.(pbipny, pbipny2, atol=1e-14))
        @test all(isapprox.(pbipnz, pbipnz2, atol=1e-14))
    end
    # 求 \bar  B_tau_xy
    bbar2 = Array{ComplexF64, 3}(undef, Nt, 2*Nx, 2*Nx)
    for bi in Base.OneTo(Nt)
        gs1 = rawgbar(allbmats1, bi-1)
        bba = bbar_from_gbar(gs1, allbmats1[bi])
        bbar2[bi, :, :] = bba
    end
    bbar = bbar_calc(ss1, allbmats1)
    @testset "B̄" begin
        @test all(isapprox.(bbar2, bbar, atol=1e-4))
    end
    # \bar nx1 = ∑_{tau,xy} \bar  B_tau_xy ∂ B_tau_xy  / ∂ nx1
    nxbar = zeros(Nx)
    nybar = zeros(Nx)#Vector{Float64}(undef, Nx)
    nzbar = zeros(Nx)#Vector{Float64}(undef, Nx)
    #for si in Base.OneTo(Nx)
    #注意虚部附加的负号
    for bi in Base.OneTo(Nt); for xi in Base.OneTo(2*Nx); for yi in Base.OneTo(2*Nx)
        nzbar += pbrpnz[bi, xi, yi, :] * real(bbar[bi, xi, yi])
        nzbar -= pbipnz[bi, xi, yi, :] * imag(bbar[bi, xi, yi])
        nybar += pbrpny[bi, xi, yi, :] * real(bbar[bi, xi, yi])
        nybar -= pbipny[bi, xi, yi, :] * imag(bbar[bi, xi, yi])
        nxbar += pbrpnx[bi, xi, yi, :] * real(bbar[bi, xi, yi])
        nxbar -= pbipnx[bi, xi, yi, :] * imag(bbar[bi, xi, yi])
    end; end; end
    #end
    #println(nxbar, nybar, nzbar)
    tapemap = pbpn_calc_tapemap(sp, nx1, ny1, nz1)
    nxbar2, nybar2, nzbar2 = meas_grad(ss1, allbmats1, sp, hscfg, nx1, ny1, nz1; tapemap=tapemap)
    @testset "n̄" begin
        @test all(isapprox.(nxbar, nxbar2, atol=1e-13))
        @test all(isapprox.(nybar, nybar2, atol=1e-13))
        @test all(isapprox.(nzbar, nzbar2, atol=1e-13))
    end
    #求sign
    sig1 = real(det(I(2*Nx) + prod(reverse(allbmats1))))
    println(sig1)
    println(-log(sig1))
    sig2 = real(det(I(2*Nx) + prod(reverse(allbmats2))))
    println(sig2)
    println(-log(sig2))
    println(-log(sig2) + log(sig1))
    println(nybar)
    #求∂ny/∂θ和∂nz/∂θ
    pnxpt = sin.(psi1) .* cos.(th1)
    pnypt = cos.(psi1) .* cos.(th1)
    pnzpt = -sin.(th1)
    pnxpp = cos.(psi1) .* sin.(th1)
    pnypp = -sin.(psi1) .* sin.(th1)
    pnzpp = zeros(Nx)
    #求∂L/∂θ
    thetabar = @. nxbar * pnxpt + nybar * pnypt +  nzbar * pnzpt
    println(thetabar)
    psibar = @. nxbar * pnxpp + nybar * pnypp +  nzbar * pnzpp
    println(psibar)
end

